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Abstract 

In this paper we present a decomposition algorithm for computation 
of the spatial-temporal optical flow of a dynamic image sequence. We 
consider several applications, such as the extraction of temporal mo¬ 
tion features and motion detection in dynamic sequences under varying 
illumination conditions, such as they appear for instance in psycholog¬ 
ical flickering experiments. For the numerical implementation we are 
solving an integro-differential equation by a fixed point iteration. 
For comparison purposes we use a standard time dependent optical 
flow algorithm, which in contrast to our method, constitutes in solving 
a spatial-temporal differential equation. 


1 Introduction 

Analyzing the motion in a dynamic sequence is of interest in many fields 
of applications, like human computer interaction, medical imaging, psychol¬ 
ogy, to mention but a few. In this paper we study the extraction of motion 
in dynamic sequences by means of the optical flow, which is the apparent 
motion of objects, surfaces, and edges in a dynamic visual scene caused by 
the relative motion between an observer and the scene. There have been 
proposed several computational approaches for optical flow computations in 
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the literature. In this paper we emphasize on variational methods. In this 
research area the first method is due to Horn & Schunck |15| . Like many 
alternatives and generalizations, the Horn & Schunck method calculates the 
flow from two consecutive frames. Here, we are calculating the optical 
flow from all available frames simultaneously. Spatial-temporal optical flow 
methods were previously studied by Weickert & Schnorr mm, mi, m 
and [2], to name but a few. However, in contrast to these paper we empha¬ 
size on a decomposition of the optical flow into components, instead of 
calculating the flow itself. 

Variational modeling of patterns in stationary images has been ini¬ 
tialized with the seminal book of Y. Meyer 123- In the context of total 
variation regularization, reconstructions of patterns was studied first in m ■ 
Here, we are implementing similar ideas as have been used before for varia¬ 
tional image denoising B HI 0 El El E9 m and optical flow decomposition 
m Eg Ea eh Eg. However, a conceptual difference is that we aim for ex¬ 
tracting temporal features, while, in all the cited papers the decomposition 
was for finding spatial components of the flow. We emphasize that the pro¬ 
posed method is one of very few variational optical flow algorithms in a 
space-time regime and within this class, this algorithm is the only spatial- 
temporal decomposition algorithm. 

The outline of this paper is as follows: In Section [2] we review the optical 
flow equation. In Section[3]we present analytical examples of the optical flow 
equation in case of illumination changes. In Section [4] we introduce the new 
model for spatial-temporal optical flow decomposition. We formulate it as a 
minimization problem and derive the optimality conditions for a minimizer. 
In Section [5] we derive a fixed point algorithm for numerical minimization 
of the energy functional. Finally in Section [6] and Section [7] we present 
experiments, results and a discussion of them. 

2 Registration and optical flow 

The problem of aligning dynamic sequences /(-,i), t G (0,1) can be formu¬ 
lated as the operator equation for finding a flow T of diffeomorphisms, 

T(-,i) : D -> D, Vi G (0,1), 

such that 


/(T(x, i),i) = /(if, 0), Mx G D and t G (0,1). (1) 

For the sake of simplicity of presentation we consider the time interval as 
the unit interval all along the paper. 

Differentiation of 0 with respect to t for a fixed x gives 

V/(T(x,i),i) • dt^ix^t) + dt/(T(x,i), t) =0, Vf G D and t G (0,1). 

( 2 ) 
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Switching from a Lagrangian to an Eulerian description we obtain the 

optical flow equation (OFE) on D: 

V/(x, t) • u(x, t) + dtf(x,t) = 0, Vf G ft and t G (0,1). (3) 

Although the derivation is based on a constant brightness assumption along 
characteristics, mathematically, ^ even makes sense under varying illumi¬ 
nation conditions. However, as we show in simple examples below, standard 
regularity assumptions on the optical flow are violated when the character¬ 
istics degenerate (collapse or originate) or when the illumination changes. 
Instead of solving ^ usually the relaxed problem is considered, which con¬ 
sists in minimizing the functional 

argmin S(u) := / (V/(x, t) • u(x, t) + dtf(x, t)) 2 dx , V£g(0, 1), (4) 

Jn 

subject to appropriate constraints. 

3 The optical flow equation in case of illumination 
changes 

In this section we are providing simple motivating examples explaining prop¬ 
erties of the solution of the optical flow equation ^ under changing illu¬ 
mination conditions. Here we are restricting attention to a spatial domain 

0 ^ Q = ( 0 , 1 ) c R. 

The following two examples simulate a day to night illumination situa¬ 
tion. The optical flow is calculated analytically, and the level lines of / are 
visualized, which represent the trajectories 4/ of constant brightness. As we 
show, smoothness of the optical flow is affected by changing illumination 
and in the first example also by joining of characteristic curves. 

Example 1. In this example the flow is not even an element of the Bochner- 
space L 2 ((0,1); i7 _1 (D)), meaning that the anti-derivative of u with respect 
to time is not in L 2 ((0,1); L 2 (D)). Because L 2 ((0,1); i7 _1 (D)) is a strict 
superset of L 2 ((0,1); L 2 (fi)), the elements have in general less regularity 
(smoothness). The next Example [2] provides a flow where / models again 
changing illumination. Here the flow is in L 2 ((0,1); i7 -1 (fi)) but not in 
L 2 ((0,1 );L 2 (fi)). We conjecture from the difference of the two examples 
that the difference in smoothness is caused by the fact that in the first 
example two characteristics are joining during the evolution. 

We consider the ID optical flow equation, to solve for u in 

d x f(x, t)u(x , t) + d t f(x, t ) = 0, V(x, t) G (0,1) x (0,1) (5) 

for the specific data 

f(x,t) = f(x)g(t), y(x, t) € (0,1) x (0,1). (6) 
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/ represents a static scene /, which is affected by brightness variations over 
time, described by g. We are more specific and take: 

f(x) = x(l — x) and g(t) = 1 — £, V(x, t) G (0,1) x (0,1). (7) 

The function / and the level lines are plotted in Figure [l] and the optical 
flow can be calculated explicitly: 

u(x,t) = V(x, t) € (0,1) x (0,1), 

indicates a transport of intensities from outside to the center 1/2. We ob¬ 
serve that u(l/2,t) and u(x, 1) are singularities of the optical flow. From 



Figure 1: f(x,t) = x(l — x)(l—t) from 0. Level lines of / are parametrized 
by 


the definition of u it follows that 

rt 


u(x,t) := j u(x,r)dT = ^ log(l ~t), V(x,t) € (0,1) x (0,1) 


and thus 


‘lli 2 ((0,l) 2 ) - f Q J I 


x 2 (l — x) z 


dx 


(1 - 2x) 2 
or in other words u L 2 ((0, l) 2 




x 2 (l — x) 2 
(1 - 2x) 2 


dx = oo, 


Example 2. This example is similar to Example [lj and simulates again a 
day to night illumination, with the difference that characteristics of / never 
join. We consider input data / of the form ^ with 

f(x) — x(l — x) and g{t) — exp | —— (1 — with some 0 < /3 < 1 (8) 

for (x,t) G Cl := (0,1/4) x (0,1). The optical flow is given by 
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Integrating this function over time gives 



u(x, r) dr 


x(l — x) 1 
l-2x /3 


(d-tf-i).. 


and consequently with 
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Figure 2: g(t) = exp |-^(1 - t )&} 


C := 


-L 


V 4 x\l - x ) 2 


dx < oo, 


we get 


Ml Ihci) = c f o ^ 2dt = 


(1 - 2x) 2 

_JC^ if/3>y 


oo else 




/3 2 V 2/3 + 1 /3 + 1 


+ 1 if p > 0 . 


This shows that u ^ L 2 (Ct) for every /? G (0, ^], but G L 2 (Ct) for all 

P G (0,1). 

The bottom line of these examples is that illumination changes, such 
as flickering, may result in singularities of the optical flow and a viola¬ 
tion of standard smoothness assumptions of the optical flow (such as u G 
L 2 ((0,1); H s {yt)) for some s > 0). The potential appearance of the singu¬ 
larities motivates us to consider regularization terms for optical flow com¬ 
putations, which allow for singularities over time, such as negative Sobolev- 
norms. 
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4 Optical flow decomposition: basic setup and for¬ 
malism 


In this paper we derive an optical flow model which allows for decomposing 
the flow field into spatial and temporal components. We consider each frame 
of the movie {/(•,£) : t G (0,1)} defined on the two-dimensional spatial do¬ 
main Q = (0, l) 2 . 

We assume that the optical flow field is a compound of two flow field 
components 


u(x, t ) =u^ (x, t) + (, x , t) = 

Vx G Q and t G (0,1). 


u { i\x,t) \ ( uf\x,t ) \ 

u£\x,t) ) ( 4 2 \x,t) J ’ 


Because there appears a series of indices and variables we specify the nota¬ 
tion in Table HJ 


The OFE-equation ^ contains four unknown (real valued) functions u^\ 
i,j = 1,2, and thus is highly under-determined. To overcome the under- 
determinacy, the problem is formulated as a constrained optimization prob¬ 
lem, to determine, for some fixed a > 0, 

argmin (ft (1) (u (1) ) + an^\u^ (9) 

subject to ©• 


Instead of solving the constrained optimization problem, we use a soft vari¬ 
ant and minimize the unconstrained regularization functional: 




2=1 


£{u^\ u*' 2 )) := /(V/ • (u (1) + u ■ 2 - ) ) + dtf) 2 dxdt with a = 


a 


( 2 ) 


( 10 ) 


a 


(i)' 


Qx (0,1) 


For the sake of simplicity of presentation, we omit arguments of the functions 
Uj ; and /, whenever it simplifies the formulas without causing misinterpre¬ 
tations. 


In the following we design the regularizes : 

• For we use a common spatial-temporal regularization functional 


for optical flow regularization (see for instance ESI): 


R<‘)(5 W):= j 


V 3 n[ 


(i) 


+ 


V 3 «? } 


dxdt, 


(ii) 


nx(o,i) 
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X = (xi,x 2 ) 

vector in two-dimensional Euclidean space 

a k = lk 

differentiation with respect to spatial variable 

II 

differentiation with respect to time 

V = (d 1 ,d 2 ) T 

gradient in space 

\7 3 = (d 1 ,d 2 ,d t ) T 

gradient in space and time 

V- = di+ d 2 

divergence in space 

v 3 - = d 1 + d 2 + d t 

divergence in space and time 

n 

outward pointing normal vector to ft 

f 

input sequence 

/(•>*) 

movie frame 


optical flow module, i = 1,2 

u — 

optical flow 

(i) 

u)> 

j -th optical flow component of the z-th mod¬ 
ule 

u(-,t) = JqU(-,t) dr 

primitive of u 

dr 

2nd primitive of u - note that dtu(-,t) =u(-,t) 


Table 1: Continuous notation. 


where v : [0, oo) —>> [0, oo) is a monotonically increasing, differentiable 
function satisfying that r —> v{r 2 ) is convex. 

According to [29] we use 


u{r) 


er + (1 — e)A 2 



Vr G [0, oo) , 


( 12 ) 


with some fixed 0 < e « 1 and A > 0. Note, that in [23 the term — 1 
is not used. However, since it is a constant, it does not influence the 
optimization. Using this term guarantees that z/(0) = 0. Moreover, 
we denote by v f the derivative of v with respect to r. 

• K® is designed to penalize for variations of the second component 
in time. Motivated by Y. Meyer’s book m, we introduce a regular¬ 
ization term, which is non-local in time. We have seen in Examples 
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[l] and [ 2 ] that u may violate L 2 -smoothness in case of changing illu¬ 
mination conditions. Variations of Meyer’s G-norm where used in 
energy functionals for calculating spatial decompositions of the opti¬ 
cal flow mm- It is a challenge to compute the G-norm efficiently, 
and therefore workarounds have been proposed. For instance m pro¬ 
posed as an alternative to the G-norm the following realization for the 
H~ x - norm: For a generalized function u : (0,1) —)> R, they defined 
IMlif- 1 — /q u^d^ufydt. Here, we use the following temporal 
H ~ x -norm for regularization: 


vJ y 2 \u^ } ) := J J u^ 2 \x,r)dr 


2 

dxdt 


Ox (0,1) 


V / (V) dMt ( l3 ) 

J=1 Ox(0,l) 


To see the analogy with the IMIjf-i-norm from |24| we note that the 
second primitive of the optical flow component vP\ as defined in Table 
0, satisfies 

Uj \x, 1) = 0, dtUj \x, 0) u^\x, 0) 0, Mj = 1,2 and x G O. 

(14) 

Then, by integration by parts it follows that 





dt 


and therefore 



(15) 


Existence of a minimizer of J 7 , defined in ^ , in an infinite dimensional func¬ 
tional analytical setting is a complicated issue. However, for some surrogate 
model, we can guarantee well-posedness by using the following lemma from 

H3: 

Lemma 1 ([!]). Let f G G 1 (fl x [0,1]). We consider t G [0,1] fixed, and 
for f := f(-,t ) we define A 0 := V/(V/) T . Then {w\,w 2 ) Ao = J Q wfA 0 w 2 
is an inner product and \w\\ 0 '■= {w,w)a 0 = Jq w t Aqw is a seminorm on 


P ■ = jX-ftfx, -ftfy? 


Let 


(16) 









1/2 i 

Then T (J = where the root of a matrix is defined via spectral 

decomposition. Moreover, 


, 1/2 - 

A) w-p 


2 

L 2 (n ; ]R 2 ) 


l|V/-«T+/t||i 2 (n) , Vn?€L 2 (f2;R 2 ). 


(17) 


Let 

-4 := ((Aq) t Aq + eld) 5, (18) 

where, Id E R 2x2 denotes the identity matrix and e > 0 is a small regu¬ 
larization parameter. A is uniformly positive definite on Q. We note that 
in comparison with p], we are not using a smoothed version of A , because 
we already made the assumption that A E C°(fl x [0,1]; R 2x2 ) (because 
/ E C 1 (fl x [0,1])). Using the notation 


$=A 2 p, 


(19) 


we assume 


l|V/-^ + /t|| L2(n) 


w — 


A 


( 20 ) 


After choosing a fixed e > 0 we consider minimization of the surrogate 
model, consisting in minimizing the functional 


X s (u^ l \u^) := £ s (u^ 1 \u^ ) ) + q/Wt^M^M), 


2=1 


^(xx^, id 2 )) := / (id 1 ) + xx (2) ) — 0 




■■ (2) 

dxdt with 
A cd 1 ) 


( 21 ) 


Ox (0,1) 


over the Bochner-Space 

x := |f7 (1 > g W 1,2 (p, x (0,1);R 2 ) : u {1] (-,T) = o} x L 2 (fl x (0,1);R 2 ). 

We are proving that the surrogate model attains a minimizer on X. Note 
that W 1,2 (fl x (0,1); R 2 ) is the space of vector valued, weakly differentiable 
functions with respect to space and time, and also recall that L 2 ((0,1); L 2 (fl; R 2 )) 
L 2 ((0,1) x fl;R 2 ). 

Theorem 1. X s attains a minimizer on X. 


Proof. The proof is done by several estimates: 

• Because (0, 0) E X it follows that C := 7^(0,0) = 



dt < oo. 
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Now, let 


(Un\u'n > ) 


be a minimizing sequence of then 



2 

A 


d t<C, 


d t < -n^iAA) < 


c 

eaW’ 



< 

L 2 (Ox(0,l);R 2 ) 


c 


• For the surrogate model ||-||^ is a norm, which is equivalent to the 

L 2 (fl; R 2 )-norm, and thus (un^+uffl) is uniformly bounded in L 2 (( 0,1); L 2 (fl)), 
and therefore admits a weakly convergent subsequence with limit </> G 
L 2 ((0,1); L 2 (fl)). This subsequence, in turn has a subsequence (which 
is also denoted by - n , for the sake of simplicity of notation) such that 

also (dfUn )) and ( u n are weakly convergent in L 2 (fl x (0,1); R 2 ), to 
some jl and respectively. 

• Let ( G L 2 (Q x (0,1); R 2 ) be arbitrary, then (x, t) -G ((x, r)dr G 
L 2 (fl x (0,1);R 2 ). Using the weak convergence of the subsequence it 
follows that 


[ • ( dxdt = f 0(x, t) • / (f(x, r)drdxdt 

Jox(0,l) Jox(0,l) Jt 

— lim [ (uffl + uff)(x, t) • [ ((x,r)drdxdt 
n ^°°Jnx( o,i) A 

[ [ (^4^ + ^)dr • (f(x,t)dtdx 

VOx(0,1) Jo 


= lim 

n—>• oo 


,/ 

5 JO: 


= lim 

7Wo ° /Ox(0,1) 


This means that (iz^) converges weakly to (f> — and in particular 

^ is uniformly bounded, let us say by D. The last item 

L 2 (Ox (0,1)) 

shows that that (L?^ 1 ^) converges to 0 — in a distributional sense. 


• C dtdx + / 
Jo 


Ox (0,1) 


i/j • £ dtdx 


( 22 ) 


10 



From integration by parts it follows that if Un\^T) = 0 that 


u, 


<i) 


/Ox (0,1) 

f d t u^\x,t) • f f u^\x,r)dr\ dtdx 
Jnx( 0 , 1 ) \J o J 


dxdt 


Jj 

d t Un ) (x, t ) 

2 dtdx [ 

/ Un\x, r) dr 

dtdx < D\l■ 

V Jnx( o,i) 


J Qx(0,l) 

Jo 

V 


C 


ea 


(i) 


As a consequence we have (j) — dt'ip G L 2 (fl x (0,1); R 2 ). 

iws from the lower semicontinr 
dt and 7 Z^(u^) for i = 1, 2. 


(23) 


• The final results follows from the lower semicontinuity of the function- 

i ^ _ 2 

als u f Q u — < 

□ 

We emphasize that the critical issue in the above proof is the coercivity, 
which could be enforced by various other modifications of the functional J 7 
as well. Another possibility, to the one mentioned above (which we feel to 
be less elegant) is to add a small regularization term for the L 2 -norm of 
to the functional J -. Since the use of surrogate models has only marginal 
impact on the numerical reconstructions we ignore their use in the following. 

4.1 Energy functional and minimization 

We are determining the optimality conditions for minimizers of J 7 introduced 


in (10). Necessary conditions for a minimizer are that the directional deriva¬ 


tives vanish for all 2-dimensional vector-valued functions h^ x (0,1) 

R 2 , j — 1,2. To formulate these conditions we use the simplifying notation: 

(£, J 7 ) := (5, T)(u^M 2) ), n® := K^{u^) and res = V/-(?2 (1) +tJ 2) )+d t f 

Therefore the directional derivative of T in direction at u = (v^K u l2 ' > ) 
is given by: 

(d nu) T)u = hrn- 1 --- 1 -= 0 

where 5ij — 1 for i — j and zero else is the Kronecker symbol. The gradient of 


the functional J 7 from (10) can be determined by calculating the directional 
derivatives of £ and separately. 


• The directional derivative of £ in direction h^ is given by 

/ res(V/ • hJ^) dxdt. 

fix (0,1) 


(24) 
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• The directional derivative of 'RS 1 ') at in direction /v 1 * is determined 
as follows: Let us abbreviate for simplicity of presentation 


v := v 


V 3 < 


(i) 


+ 


v 3 4 1} 


u’ := u’ 


V 3 «i 


(i) 


+ 


v 3 4 1} 


then the directional derivative of in direction /T 1 ) (which we as¬ 
sume to have compact support in x (0,1)) at u ^ is given by 


v ^ (1) (^ (1) + shW) -nW 

ftn * = >™ 0 - 1 - 7 - L - 


= - 2 


[ V 3 • (Vv 3 r4 1)N ) + V 3 • (Vv 3 4 1} ) h^dxdt, 

Jnx( 0,1) V J V 7 


(25) 


where integration by parts is used to prove the final identity. 

• The directional derivative of is derived as follows: 

i-c y “'\u y “' + sn y “’ i - 'K y “> n . 2 
= lim-^-= 2 


(^ ( 2 )7e (2) )n (2) = lim 


2 r 

T-l 


=1 JQx( 0,1) 


dxdt. 


(26) 


Moreover, it follows by integration by parts of the last line of (26) with 
respect to t that 


(8 SB ,K< 2| )fi< 2) = -2]T f 

j=1 Jnx( o,i) 


u[ \^dxdt. 


J 3 


(27) 


Now, because of (25) and (24) it follows that the minimizer u^\ i — 1, 2 


has to satisfy for every j = 1,2, 

djfiyf • (S (1) + U {[2) ) + dtf ) - « (1) v 3 • (VV 3 ^ 1} ) = 0 in n x (0,1), 

9ft u ^ — 0 in dfl x (0,1), 


— 0 in x {0,1} . 


(28) 

->( 2 ) 

Since equations (24) and (27) hold for all h) , it follows that for every 

J = 1,2, 


djfiyf • (u^ + i^ 2 )) + dtf ) — )= 0 in x (0,1) 


Thus the optimality conditions for a minimizer are (28) and 


(29) 


The 


solution of Equation (28) has to be understood in a weak sense. The opti¬ 


mality condition is formal and would be exact if T would be coercive. As 
we demonstrated above there are several possibilities of surrogate models, 
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/ = £ r M x ^ xT 

input sequence 

= u^(r, s, t; k) € Wi MxNxTxKx2 

discrete optical flow approximating the 
continuous flow a t (yfir[, aTT’ T^l) 

4 

finite difference approximation in direction x & 

$ 

finite difference approximation in direction t 

— M-\ ’ — N- 1 an< ^ — T—l 

Discretization 

uf\r,s,t-,k ) = A* Xr=i u f ) ( r > s > T ; * 0 > 3 = 1,2 

finite difference approximation of 

(r, s, t ; fc) = - At Xr=t “f \ r > 

finite difference approximation of S(-,t) 


Table 2: Discrete Notation. 

which provide coercivity. The surrogate model outlined above leading to 
this optimality condition would require a Dirichlet boundary condition for 
u W at D x { 1 }. If we would indeed complement T by a 3 H'uW 
with small 0 ^ 3 , then the boundary conditions would in fact be the natural 
ones. 


L 2 (Ox (0,1)) ’ 


5 Numerics 

In this section we discuss the numerical minimization of the energy func¬ 


tional T defined in (10). Our approach is based on solving the optimality 
conditions for the minimizer u^\u^ from (28), (29) with a fixed point it¬ 


eration. We call the iterates of the fixed point iteration u^\x,t\k), where 
k = 1, 2,..., K denotes the iteration number and K is the maximal num¬ 
ber of iterations. We summarize all the iterates of the components of flow 

(i) 

functions u K - in a tensor of size M x N x T x K . In this section we use the 
notation as reported in Table [2j 

For every tensor H — (H(r, s,t)) E ]& MxNxT (representing all frames of 
a complete movie) we define the discrete gradient 


V%H(r,s,t) = (dtH(r,s,t),d%H(r,s,t),d?H(r,s,t)) T 

for (r, s, t) € M} x {1,...,N} x ,T}, 
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where 


^ _H(r + l,s,t) - H(r-l,s,t) 


d\H(r,s, t) 


2 A* 


if 1 < r < M, 


a?g(r, + lf , < , < N , (30) 


2 A y 

^hu, ,x _H(r,s,t + 1) - H(r,s,t - 1) 


d?H(r,s,t) 


2A t 


if 1 < t < T. 


Let us notice that r,s,t are discrete indices corresponding to the points in 
space r A x , sA y and in time tA t of a continuous function H. Moreover, we 
define the discrete divergence, which is the adjoint of the discrete gradient: 
Let (Hi, H 2 , H 3 ) T (r, s,t), then 

V 3 • (Hi,H 2 , H 3 ) t = d\H, + d\H 2 + d\ l H 3 . (31) 

The realization of the fixed point iteration for solving the discretized equa¬ 


tions (28) and (29) reads as follows: 


• Initialization for A: = 0: we initialize to 0 ^^(O), i/( 2 )(0) E ]^ MxNxKx2 1 
Note, we leave out all indices denoting space and time positions. 

• k -> k + 1 and k < K: let v'(k) := ^ , (|V^[ 1 ) (fc )| 2 + |V^ 2 ) (fc)| 2 ), 
then 


(k + 1 ) — u\ LI (k) 

~ a; 


a), 


= \/ h 3 -(u / (k)V h 3 u\ L \k)) 
d h if 


a 


(i) 


dif (u^\k + 1 ) + u[ 

+ 9 2 f («2 \k) + u 2 \k)^J + dtf 


(32) 


^(fc + l) v£>(k) = . { ^ m h u (i) m 


,(i) 


A 


dif 


a 


(i) 


dif («i 1 ) (A: + 1 ) +u[ 

+ d 2 f («2 \k + 1) + u. 


+ dif 


u^\k + 1 ) — uY\k ) 

“ a: 


( 2 )/ 


(33) 


= uf\k) 

dif 


a 


( 2 ) 


[dif (u^/c + 1) + ui\k + 1)) 
-\~d 2 f («2 ^ (k + 1 ) + ui. 


+m 


( 34 ) 
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and 


u ?\ k + !) - 4 2) ( fc ) = 


A t 


d'lf 


( 2 ) 


dif (u^\k + 1) + u{\k + 1 

+^ 2 / + 1) + + 1)) + ®tf 

(35)' 

where A r is the step size iteration parameter, which has been set 
to 10 -4 . In order to improve stability of the algorithm and ensure 
convergence we use a semi implicit iteration scheme as proposed in 


Indeed, let us notice that in (32), (33),(34), (35) one of the 
components of the two flow fields Hy\ on the right hand side 


is evaluated for iteration index k + 1. The system (32), (33), (34), 


(35) can be solved efficiently using the special structure of the matrix 
equation (see |2Hl ESI). The matrix equation then is a sparse block 


tridiagonal matrix.This type of matrices allows for efficient estimation 
of a solution through decomposition and parallelization methods 


The iterations are stopped when the Euclidean norm of the relative error 
\uf(-,k)-uf(-,k + l)\ 

drops below the precision tolerance value tol = 0.05 for at least one of the 
component of u W and one of u ^. The typical number of iterations is much 
below 100. 


6 Experiments 

In this section we present numerical experiments to demonstrate the po¬ 
tential of the proposed optical flow decomposition model. In the first two 
experiments we use for visualization of the computed flow fields the standard 
flow color coding [8]. The flow vectors are represented in color space via the 
color wheel illustrated in Figure [3| For the third and fourth experiment we 
use a black and white visualization technique: Black means that there is 
no flow present and the gray-shade is proportional to the flow magnitude. 
In order to compare the optical flow computations quantitatively the inten¬ 
sity values of / have been scaled to the range [0,1]. The used parameters 
are reported for each experiment separately, except for the discretization 
parameters A x ,A y , At, defined in Section [5j In this work we consider the 
following four dynamic image sequences: 
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Figure 3: Color Wheel. 

• The first experiment is performed with the video sequence from m 
(available at http://of-eval.sourceforge.net/) which consists of 
forty-six frames showing a rotating sphere with some overlaid patterns. 
The analytical results from Appendix [A] in ID suggest that the inten¬ 
sity of the u component increases monotonically with increasing 
rotational frequency over time. We verify this hypothesis numerically, 
now in higher dimensions. We simulate in particular two, four and 
eight times the original motion frequency;: In order to do so, we du¬ 
plicate the sequence periodically, however consider it to be recorded 
in the same time interval (0,1). The flow visualized in Figure [6] is 
the one between the 16th and the 17th frame of every sequence. We 
study the behavior of the sphere at different motion frequencies with 
the same parameter setting a W = 1, a^ The numerical results 

confirm the ID observation that for high rotational movement u is 
dominant (cf. Figures ^ and is always 20% of u in particular 
u W and u^ cannot be completely separated. 




•v 


Figure 4: u^ at different frequencies of rotations: 2, 4 and 8x faster than 
the original motion frequency, a^ 1 ) = 1, a^ The intensity of u^ 

increases when the frequency of rotation is increased. 


> *•• ** x \ 

</• .* *V 

v v- ••J 


• The second experiment concerns the decomposition of the motion in a 
dynamic image sequence showing a projection of a cube moving over an 
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oscillating background. The movie consists of sixty frames and can be 
viewed on the web-page http: //www. esc. univie. ac. at/index. php? 
page=visualattention, 

The background is oscillating in diagonal direction, from the bottom 
left to the top right, with a periodicity of four frames. In each frame 
the oscillation has a rate of 5% of the frame size. The flow visualized 
in Figure [6] is the one between the 20th and the 21st frame of the 
sequence. 

Applying the proposed method with a parameter setting a ^ = 10 3 , 
q/( 2 ) = 10 3 , A r = 10 -5 , and precision tolerance tol = 0.001, we notice 
that the background movement appears almost solely in u and the 
global movement of the cube appears in In Figure [6] we represent 
only flow vectors with magnitude larger than 0.05 and omit in u^ the 
part in common with u ^ for better visibility. 



Figure 5: The dynamic sequence consists of the smooth (translation like) 
motion of a cube and an oscillating background. The oscillation has a pe¬ 
riodicity of four frames and takes place along the diagonal direction from 
the bottom left to the top right, moving at a rate of 5% of the frame size 
in each frame. The proposed model decomposes the motion, obtaining the 
global movement of the cube in u ^ (left) and the background movement 
exclusively in u^ (right). 


• In the third experiment the original movie consists of thirty-two frames 
and can be viewed together with the decomposition result on the web¬ 
page http : //www.esc.univie.ac.at/index.php?page=visualattention, 
The flow is decomposed into two components. The first part shows the 
movement of a Ferris wheel and people walking. The second part shows 
blinking lights and the reflection of the wheel. The flow visualized in 
Figure [6] is the one between the 4th and the 5th frame of the sequence 
with a parameter setting a W = 1, a In order to improve the 
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Figure 6: u^: Movement of a Ferris wheel and people walking in the fore¬ 
ground (top left), u consists of blinking lights and the reflections of the 
wheel (top right). The third image (bottom) is a reference frame. 


visibility we represent only flow vectors with magnitude larger than 
0.18 and we omit for u the part in common with 

• The fourth example is flickering. In a standard flickering experi¬ 
ment, the difference in human attention is investigated by inclusion 
of blank images in a repetitive image sequence. Although, in gen¬ 
eral, these blank images are not deliberately recognized, they change 
the awareness of the test persons. J. K. O’Regan [22] states that 
“Change blindness is a phenomenon in which a very large 
change in a picture will not be seen by a viewer, if the 
change is accompanied by a visual disturbance that prevents 
attention from going to the change location”. They test data 
from http://nivea.psycho.univ-paris5.fr, was used for our sim¬ 
ulations. The proposed optical flow decomposition is able to detect 
regions, which also humans can recognize, but standard optical flow 
algorithms fail to: To show this, the input sequence is composed by 
four frames consisting of Frame 1, a blank image, Frame 2 and again 
an identical blank image (see Figure [8] (top)). This sequence is then 
aligned periodically to a movie. We interpret the movie as a linear in¬ 
terpolation between the frames. We test and compare Horn-Schunck, 
Weickert-Schnorr and the proposed algorithm. 

We set the smoothness parameter a ^ to a value of one for all the 
methods. Moreover, for our approach we set a = 1. For Horn- 
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Figure 7: Result with Horn-Schunck. 


Schunck we visualize the flow field in Figure [7| This flow is the one 
between the blank frame and the slightly changed frame, which exceeds 
a threshold of 3.9. The results obtained by applying Weickert-Schnorr 
and the field of our approach, respectively, are small in magni¬ 
tude. Therefore, we do not visualize them. This behavior is coherent 
with the motivation of the Weickert-Schnorr method to produce an 
optical flow that is less sensitive to variations over space and time. 
Finally, we visualize in Figure [5] (down right) the u^ flow field for 
the proposed approach. For the visualization we omit all vectors with 
magnitude lower than 0.18. In order to make transparent the result, 
we show in Figure [8] (down left) the difference between the two frames 
of the sequence containing information (see Figure [8] (top)). In this 
experiment, we notice that the u W component is negligible, instead 
u detects the areas affected by change of intensities (see Figure ^ 
(down right)). 

Additional Information 

In the following, we show the capacity of our model to extract different 
information compared to standard optical flow algorithms. The current lit¬ 
erature focuses on average angular and endpoint error [8] in order to compare 
optical flow algorithms. Our model extracts information, that is neglected 
by standard algorithms. Such difference can be shown through a quanti¬ 
tative comparison of models. For this purpose, we use well-known test se¬ 
quences from the Middlebury database http://vision.middlebury.edu/ 
flow/, and evaluate the residual of the optical flow constraint. 

In order to understand how much information our method is capable to 
extract from an entire dynamic sequence, we calculate the residual squared 
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Figure 8: The two frames of the flickering sequence containing information 
(top), the difference between these two frames (down left), and the u^ flow 
field resulting from the proposed approach (down right). As predicted in 
Section [ 3 ] and Appendix [a] the u^ component is negligible, instead 
detects the change of intensity across the blank sheet. 

over space and time: £(u^\u^) as in ( [T()| ) and compare it with the squared 
residual over space and time of the Weickert-Schnorr method [28, 29]. We 
use the parameter settings a ^ = 100 (a = in Weickert-Schnorr) and 
a tolerance tol = 0.01, in order to have a good comparison of the two 
methods. Again the residual is smaller for the proposed method as shown 
in Table HI 

The smaller value of the residual are due to the fact that we are calcu¬ 
lating a minimizer from a larger space of optical flow components. However, 
with this approach we cannot observe oversmoothing of the flow. 



Weickert-Schnorr 

Proposed model 

Hamburg Taxi 

1374.9 

1021 

Rubber Whale 

4459.7 

3046.8 

Hydrangea 

8533.3 

7647.2 

DogDance 

9995.4 

8217.6 

Walking 

8077.5 

5944.3 


Table 3: Comparison of squared residuals over space and time £ between 
Weickert-Schnorr and the proposed method. 
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7 Conclusion 


We present a new optical flow model for decomposing the flow in spatial 
and temporal components of different scales. A main ingredient of our work 
is a new variational formulation of the optical flow equation. Finally, appli¬ 
cations (some of them from psychological experiments) are considered and 
analyzed analytically and computationally. 
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A Optical flow decomposition in ID 


In order to make transparent the features of our decomposition we study 
exemplary the ID case again. From regularization theory (see e.g. [23]) 
we know that the minimizers (uQ\ ix^), for a = 0, are 

converging to a solution of the optical flow equation which minimizes 

( 2 ) 

1Z = for a = lim — j-r > 0 . 

a—>-0 cdB 

Such a solution is called 1Z minimizing solution. Note that by the ID sim¬ 
plification the modules i — 1,2, are real valued functions. 

We calculate the decomposition for the optical flow equation ([5]) , for the 
specific test data <§• The regularized solutions approximate the 1Z min¬ 
imizing solution, and thus these calculations can be viewed representative 
for the properties of the minimizer of the regularization method. For these 
particular kind of data the solution of the optical flow equation is given by : 


u(x, t) = — 


f(x) d t g(t ) _ 3t(log g)(t) 


(36) 


d x f(x) g(t ) <9* (log f)(x) 

Let us assume that (lo gg)(t) — (logg)(0) can be expanded into a Fourier 


sm-series: 


rt 00 

(log 0 )(*) - (logs)(0) = / (logsXO dr = W sin(nTrf). (37) 

Jo n=l 
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Moreover, we assume that 1 /d x (\og f)(x) can be expanded in a cos-series: 

J oo 

-~-= > fm cos(mnx). 

4 (log /)(*) ^ 

Then 

Hncr n\(i\ — rinanVO) /-.n / 0 x 

- 2(x, t ) t ) + ?r 2 '(£, t). 


(38) 


(39) 


<9* (log/) (a;) 

Inserting this identity in the regularization functional TZ(u^\ u i2> ) = J{d x v {l> ) 2 + 
(dtu 11 )) 2 + a (u (2) 'j 2 dxdf , we remove the u (2> dependence, and we get 


7^(n (1) ) := f {d xt u^) 2 +(d tt u^) 2 +ot 

J( 0,l)x(0,l) 


/ (logg)(t) - (logff)(0) 
V <9* (log/) (a;) 


\ 2 

+ j dxdt, 


where we enforce the following boundary conditions on {(I 1 ): From the defini¬ 
tion of n! 1 ) it follows that u^\x, 0) = 0. Secondly, we enforce u ( 1 - 1 (.x, 1) = 0. 
In fact, the assumption is reasonable, because when the series 9™ I s 

absolutely convergent, (log g)(t) — (logg)(0) = 0 (cf. ©), which implies 
that 1) + u( 2 \x, 1) = 0 (cf. (|37|)). 

By substituting the relation between u^ and u ^ we reduce the con¬ 
straint optimization problem to an unconstrained optimization problem for 
and the minimizer solves the partial differential equation 

dttxxU {l) + dttttu (1) +a( — (logg)(°) + ~(i)\ = o in (0,1) x (0,1), 

V 9* (log f)(x) J 

together with the boundary conditions: 

dttu ^ — u ^ — 0 on (0,1) x {0,1} , 
d x dttu (1) = 0 on {0,1} x (0,1). 


(40) 


The boundary conditions dttu^ == 0 on (0,1) x {0,1} , and d x dttu ^ = 
0 on {0,1} x (0,1) appear as natural boundary conditions, when weak so¬ 
lutions are considered. 

Now, we substitute w := dttu^\ and we get the following system of 
equations 


d xx w + d u w = -a 


^ (log 9)(t) - (logffKO) 


(log f)(x) 
w — 0 on (0,1) x {0,1} , 
d x w — 0 on {0,1} x (0,1). 


+ in (0,1) x (0,1), 


( 41 ) 


and 


rt rr r 1 rr 

u^\x,t) — / / w(x, r)drdr — t / / u/(x,f)dfdr. 

Jo Jo Jo Jo 
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w can be expanded as follows: 

oo 

w(x,t) — WmnCOs(mnx) sm(nnt), 

171 , 71 =0 

and we expand u W in an analogous manner: 

oo 

u^\x,t) = u$ n cos(m7rx) sm(ri7rt) , 

m,n=0 


such that 


VJr 


= -n 2 7r 2 n^, Vm, n G N 0 . 


Thus it follows from (41) that 


w mn (m 2 + n 2 )ir 2 = a (u^ n + f m g^j , Vm,n <E N 0 . 


(42) and (43) imply that 


u (!) = _ 

mn 


a 


a + 7 r 4 (m 2 + n 2 )n 2 


fm9n , Vm,n G Nq. 


Now, consider a specific test example g(£) = exp j | for 

G N. Then, froi 
spondingly we have 


no G N. Then, from (1361) it follows that u(x,t) — — c< **( nQ7 ^) and 

- v ' ^(log/)(a?) 7 


(42) 

(43) 

(44) 

some 

corre- 


(log g)(t) - (log g)(0) = 


sin(no7rt) 

n 0 7r 


V ^o s in(n7rt). 

“ n 0 7r 
n=l 


In this case it follows from (44) that 


u, 


(!) = -- 


a 


J nn o 


a + 7r 4 (m 2 + njj)nQ no7r 


,/m- 


In the case of flickering data /, u^ is pronounced (if no is large Umh ~ 
0) while in the quasi-static case is dominant. Moreover, we also see 
that spatial components belonging to Fourier-cos coefficients with large m 
are more pronounced in the u^ component, and the spatial and temporal 
coefficients always appear in both components. In particular this means 
that a threshold has to be set, to assign them to the first or second module. 
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